rm(list=ls())
setwd("/GhanaPaper_Governance/ReplicationFiles")

library(arm)
library(haven)

set.seed(1)

# Turning off scientific notation
options(scipen = 999)

# Data
data_full <- read.csv("/Data/pol_allline_r.csv")


#####Figure 3##########
# Estimated coefficients and standard errors come from 
# Table C1. Refer to appendix_c.do.

jpeg("Figure 3.jpeg")

coefplot(c(0.3886, 0.1923),
         c(0.1010, 0.0990),
         varnames = c("Control", "GAS"),
         col.pts = "red",
         main = "Partisan Bias", var.las = 0, xlim = c(0, 0.8), vfont = c("sans serif", "plain"), cex.pts = 1.5, cex = 1.3)
text(0.4394, 1.2, "0.4394")
text(0.5989, 2.2, "0.5989")


coefplot(c(0.4394, 0.5989),
         c(0.1167, 0.1015),
         add = T, cex.pts = 1.5)
text(0.3886, 0.9, "0.3886", col = "red")
text(0.1923, 1.9, "0.1923", col = "red")

text(0.3, 0.7, "Baseline vs.", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.45, 0.7, "Endline", col = "red", vfont = c("sans serif", "plain"), cex = 1.1)

dev.off()

#####Figure 4##########
## DiD Estimates for GAS

# Subsets of the data
endline <- data_full[which(data_full$survey_wave == 2),]
baseline <- data_full[which(data_full$survey_wave == 0),]
midline <- data_full[which(data_full$survey_wave == 1),]

con_endline <- endline[which(endline$trt == 0),]
gas_endline <- endline[which(endline$trt == 1),]

con_baseline <- baseline[which(baseline$trt == 0),]
gas_baseline <- baseline[which(baseline$trt == 1),]

con_midline <- midline[which(midline$trt == 0),]
gas_midline <- midline[which(midline$trt == 1),]

# Estimates from the List Experiment in Control Districts
dif_con_baseline <- mean(con_baseline$surveyexp3[con_baseline$list_3==1], na.rm = TRUE) - mean(con_baseline$surveyexp3[con_baseline$list_3==0], na.rm = TRUE)
dif_con_endline <- mean(con_endline$surveyexp3[con_endline$list_3==1], na.rm = TRUE) - mean(con_endline$surveyexp3[con_endline$list_3==0], na.rm = TRUE)
dif_con_midline <- mean(con_midline$surveyexp3[con_midline$list_3==1], na.rm = TRUE) - mean(con_midline$surveyexp3[con_midline$list_3==0], na.rm = TRUE)

# Estimates from the List Experiment in GAS Districts
dif_gas_baseline <- mean(gas_baseline$surveyexp3[gas_baseline$list_3==1], na.rm = TRUE) - mean(gas_baseline$surveyexp3[gas_baseline$list_3==0], na.rm = TRUE)
dif_gas_endline <- mean(gas_endline$surveyexp3[gas_endline$list_3==1], na.rm = TRUE) - mean(gas_endline$surveyexp3[gas_endline$list_3==0], na.rm = TRUE)
dif_gas_midline <- mean(gas_midline$surveyexp3[gas_midline$list_3==1], na.rm = TRUE) - mean(gas_midline$surveyexp3[gas_midline$list_3==0], na.rm = TRUE)

## DiD estimates: 
# GAS Endline - Midline
did_gas1 <- (dif_gas_midline - dif_con_midline) - (dif_gas_baseline - dif_con_baseline)
# GAS Endline - Midline 
did_gas2 <- (dif_gas_endline - dif_con_endline) - (dif_gas_midline - dif_con_midline)
# GAS Endline - baseline 
did_gas3 <- (dif_gas_endline - dif_con_endline) - (dif_gas_baseline - dif_con_baseline)

## Standard Errors for the DiD estimator (GAS)
## Bootstrap standard errors, clustering at the TRIPLET level

r <- 1000                       # number of bootstrap samples
remove <- c(0)

# Vectors to hold estimates for each bootstrap sample
dif_con_baselinev <- rep(0,r)
dif_con_endlinev <- rep(0,r)
dif_con_midlinev <- rep(0,r)
dif_gas_baselinev <- rep(0,r)
dif_gas_endlinev <- rep(0,r)
dif_gas_midlinev <- rep(0,r)

did_gasv1 <- rep(0,r)
did_gasv2 <- rep(0,r)
did_gasv3 <- rep(0,r)

ones <- rep(1,r)                # vector of ones

for (i in (1:r)){
  b <-c()
  b1 <- c()
  b2 <- c()
  b3 <- c()
  # Creating each chunk
  ### Resampling Districts per Survey wave
  district_bas <-sample(1:length(unique(baseline$triplet,incomparables = FALSE)), replace=TRUE)
  district_mid <-sample(1:length(unique(midline$triplet,incomparables = FALSE)), replace=TRUE)
  district_end <-sample(1:length(unique(endline$triplet,incomparables = FALSE)), replace=TRUE)
  ldistrict_bas <- length(district_bas)
  ldistrict_mid <- length(district_mid)
  ldistrict_end <- length(district_end)
  for (h in 1:ldistrict_bas){
    bdist <- baseline[which(baseline$triplet == district_bas[h]),]
    b1<-rbind(b1,bdist)
  }
  for (j in 1:ldistrict_mid){
    bdist2 <- midline[which(midline$triplet == district_mid[j]),]
    b2<-rbind(b2,bdist2)
  }
  for (k in 1:ldistrict_end){
    bdist3 <- endline[which(endline$triplet == district_end[k]),]
    b3<-rbind(b3,bdist3)
  }
  b <- rbind(b,b1,b2,b3)
  # Data subsets according to Survey Wave and Treatment Assignment
  endline_b <- b[which(b$survey_wave == 2),]
  baseline_b <- b[which(b$survey_wave == 0),]
  midline_b <- b[which(b$survey_wave == 1),]
  
  con_endline_b <- endline_b[which(endline_b$trt == 0),]
  gas_endline_b <- endline_b[which(endline_b$trt == 1),]
  
  con_baseline_b <- baseline_b[which(baseline_b$trt == 0),]
  gas_baseline_b <- baseline_b[which(baseline_b$trt == 1),]
  
  con_midline_b <- midline_b[which(midline_b$trt == 0),]
  gas_midline_b <- midline_b[which(midline_b$trt == 1),]
  
  # List Esperiment 3 - Differences within each survey wave and Audit Score type  
  # Control
  dif_con_baselinev[i] <- mean(con_baseline_b$surveyexp3[con_baseline_b$list_3==1], na.rm = TRUE) - mean(con_baseline_b$surveyexp3[con_baseline_b$list_3==0], na.rm = TRUE)
  
  dif_con_endlinev[i] <- mean(con_endline_b$surveyexp3[con_endline_b$list_3==1], na.rm = TRUE) - mean(con_endline_b$surveyexp3[con_endline_b$list_3==0], na.rm = TRUE)
  
  dif_con_midlinev[i] <- mean(con_midline_b$surveyexp3[con_midline_b$list_3==1], na.rm = TRUE) - mean(con_midline_b$surveyexp3[con_midline_b$list_3==0], na.rm = TRUE)
  
  # GAS
  dif_gas_baselinev[i] <- mean(gas_baseline_b$surveyexp3[gas_baseline_b$list_3==1], na.rm = TRUE) - mean(gas_baseline_b$surveyexp3[gas_baseline_b$list_3==0], na.rm = TRUE)
  
  dif_gas_endlinev[i] <- mean(gas_endline_b$surveyexp3[gas_endline_b$list_3==1], na.rm = TRUE) - mean(gas_endline_b$surveyexp3[gas_endline_b$list_3==0], na.rm = TRUE)
  
  dif_gas_midlinev[i] <- mean(gas_midline_b$surveyexp3[gas_midline_b$list_3==1], na.rm = TRUE) - mean(gas_midline_b$surveyexp3[gas_midline_b$list_3==0], na.rm = TRUE)
  
  #Difference in Difference Estimators 
  # Midline - Baseline
  did_gasv1[i] <- (dif_gas_midlinev[i] - dif_con_midlinev[i]) - (dif_gas_baselinev[i] - dif_con_baselinev[i])
  # Endline - Midline
  did_gasv2[i] <- (dif_gas_endlinev[i] - dif_con_endlinev[i]) - (dif_gas_midlinev[i] - dif_con_midlinev[i])
  
  # Endline - Baseline
  did_gasv3[i] <- (dif_gas_endlinev[i] - dif_con_endlinev[i]) - (dif_gas_baselinev[i] - dif_con_baselinev[i])
  
  # Endline - Baseline for CSO
  did_csov[i] <- (dif_cso_endlinev[i] - dif_con_endlinev[i]) - (dif_cso_baselinev[i] - dif_con_baselinev[i])
  
}
# Standard Errors: 
se_v2 <- rep(0,3)

# SE For Outcomes
# Midline - Baseline
se_v2[1] <- sqrt(sum((did_gasv1 - (ones*did_gas1))^2, na.rm =TRUE)/(r-1))
# Endline - Midline
se_v2[2] <- sqrt(sum((did_gasv2 - (ones*did_gas2))^2, na.rm =TRUE)/(r-1))
# Endline - Baseline
did_bse <- sqrt(sum((did_gasv3 - (ones*did_gas3))^2, na.rm =TRUE)/(r-1))
se_v2[3] <- did_bse

## DiD estimates and Standard Errors
dif_v <-rep(0,3)
dif_v[1] <- did_gas1
dif_v[2] <- did_gas2
dif_v[3] <- did_gas3

cbind(dif_v, se_v2)

## Generating the Plot:
jpeg("Figure 4.jpeg")

coefplot(c(did_gas3),
         c(did_bse),
         varnames = c("GAS"),
         main = "Difference-in-differences Estimates",
         ylab = "Treatment", cex.pts = 1.5, cex = 1.3, xlim= c(-0.8, 0.4), cex.var = 1.3)
text(did_gas3, 1.1, str(did_gas3))

dev.off()


#####Figure 5##########
# Estimated coefficients and standard errors come from 
# Table C1. Refer to appendix_c.do.

jpeg("Figure 5.jpeg")

coefplot(c(0.1923, 0.1989, 0.5989),
         c(0.0990, 0.1127, 0.1015),
         varnames = c("Endline", "Midline", "Baseline"),
         main = "Partisan Bias", xlim = c(-0.2, 1),
         ylab = "Rounds", cex.var = 1.3, cex.pts = 1.4, cex = 1.2, var.las = 0, col.pts = c("black", "blue", "red"))
text(0.1923, 1.1, "0.1923")
text(0.1989, 2.1, "0.1989")
text(0.5989, 3.1, "0.5989")

dev.off()


#####Figure 6##########
# Estimated coefficients and standard errors come from 
# Table C2. Refer to appendix_c.do.

jpeg("Figure 6.jpeg")

coefplot(c(0.2626, 0.1167),
         c(0.1543, 0.1391),
         varnames = c("Positive Audit Score", "Negative Audit Score"),
         main = "Partisan Bias", var.las = 0, xlim = c(-0.2, 1.2), cex.pts = 1.5, cex = 1.2, col.pts = "red")
text(0.2626, 0.95, "0.2626", col = "Red")
text(0.1167, 1.95, "0.1167", col = "Red")

coefplot(c(0.1222, 0.2959),
         c(0.1695, 0.1757), add = T,
         col.pts = "blue", cex.pts = 1.5)
text(0.1222, 1.18, "0.1222", col = "Blue")
text(0.2959, 2.18, "0.2959", col = "Blue")

coefplot(c(0.3487, 0.8463),
         c(0.1676, 0.1346),
         add = T,
         offset = 0.23, cex.pts = 1.5)
text(0.3487, 1.32, "0.3487")
text(0.8463, 2.32, "0.8463")

text(0.3, 0.5, "Baseline vs.", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.6, 0.5, "Midline vs.", col = "blue", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.85, 0.5, "Endline", col = "red", vfont = c("sans serif", "plain"), cex = 1.1)

dev.off()


#####Figure 7##########
# Estimated coefficients and standard errors come from 
# Table C3. Refer to appendix_c.do.

jpeg("Figure 7.jpeg")

coefplot(c(0.1482, 0.2194),
         c(0.1786, 0.1221),
         varnames = c("NPP", "NDC"),
         main = "Partisan Bias", var.las = 0, xlim = c(-0.4, 1.0), cex.pts = 1.5, cex = 1.2, col.pts = "red")
text(0.1482, 0.95, "0.1482", col = "Red")
text(0.2194, 1.95, "0.2194", col = "Red")

coefplot(c(-0.1232, 0.3108),
         c(0.1668, 0.1524), add = T,
         col.pts = "blue", cex.pts = 1.5)
text(-0.1232, 1.18, "-0.1232", col = "Blue")
text(0.3108, 2.18, "0.3108", col = "Blue")

coefplot(c(0.5833, 0.6065),
         c(0.2283, 0.1328),
         add = T,
         offset = 0.23, cex.pts = 1.5)
text(0.5833, 1.32, "0.5833")
text(0.6065, 2.32, "0.6065")

text(0.3, 0.5, "Baseline vs.", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.6, 0.5, "Midline vs.", col = "blue", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.85, 0.5, "Endline", col = "red", vfont = c("sans serif", "plain"), cex = 1.1)

dev.off()


#####Figure 8##########
# Estimated coefficients and standard errors come from 
# Table C4. Refer to appendix_c.do.

jpeg("Figure 8.jpeg")

coefplot(c(0.1808, 0.2205),
         c(0.1251, 0.1767),
         varnames = c("Low", "High"),
         main = "Partisan Bias", var.las = 0, xlim = c(-0.2, 1.0), cex.pts = 1.5, cex = 1.2, col.pts = "red")
text(0.1808, 0.95, "0.1808", col = "Red")
text(0.2205, 1.95, "0.2205", col = "Red")

coefplot(c(0.1488, 0.2915),
         c(0.1607, 0.1774), add = T,
         col.pts = "blue", cex.pts = 1.5)
text(0.1488, 1.18, "0.1488", col = "Blue")
text(0.2915, 2.18, "0.2915", col = "Blue")

coefplot(c(0.5324, 0.6971),
         c(0.1685, 0.1406),
         add = T,
         offset = 0.23, cex.pts = 1.5)
text(0.5324, 1.32, "0.5324")
text(0.6971, 2.32, "0.6971")

text(0.3, 0.5, "Baseline vs.", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.55, 0.5, "Midline vs.", col = "blue", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.75, 0.5, "Endline", col = "red", vfont = c("sans serif", "plain"), cex = 1.1)

dev.off()

#####Figure F1##########
# Estimated coefficients and standard errors come from 
# Table C1. Refer to appendix_c.do.
jpeg("Figure F1.jpeg")

coefplot(c(0.3886, 0.1923, 0.3545),
         c(0.1010, 0.0990, 0.0933),
         varnames = c("Control", "GAS", "CSO"),
         col.pts = "red",
         main = "Partisan Bias", var.las = 0, xlim = c(0, 0.8), vfont = c("sans serif", "plain"), cex.pts = 1.5, cex = 1.3)
text(0.4394, 1.2, "0.4394")
text(0.5989, 2.2, "0.5989")
text(0.4470, 3.2, "0.4470")

coefplot(c(0.4394, 0.5989, 0.4470),
         c(0.1167, 0.1015, 0.1095),
         add = T, cex.pts = 1.5)
text(0.3886, 0.9, "0.3886", col = "red")
text(0.1923, 1.9, "0.1923", col = "red")
text(0.3545, 2.9, "0.3545", col = "red")

text(0.3, 0.7, "Baseline vs.", vfont = c("sans serif", "plain"), cex = 1.1)
text(0.45, 0.7, "Endline", col = "red", vfont = c("sans serif", "plain"), cex = 1.1)

dev.off()
